% The coefficient array obtainment for N dimensional Chebyshev Polynomials Gauss by Chebyshev Integration and Smolyak sparse grid technique
% It is noted that this m file need to be modified manually!
% a: the lower limit
% b: the upper limit
% Polynomial_Order: the total degree of N dimensional Chebyshev Polynomials
% func: the function handle
function [Coefficient_Set,Degree_Nd_Set] = Surrogate_Model_Nd(a,b,Polynomial_Order,func)

	Ndim = 6;																	% the number of dimensions
	N_kLevel = 10;																% the number of k - Level
	Num_monomials = nchoosek(Polynomial_Order+Ndim,Ndim);                       % the number of monomials
	NumRange_Gpoints_1d = 1 : Ndim+N_kLevel;                                    % the range of 1D Gauss integration point amount,this range cannot contain the element '0'
	DegRange_1d = 1 : Polynomial_Order+1;                                        % the degree range of a variable, such as x, y, z
	[Num_Sequence_No1d, Num_Sequence_No2d,Num_Sequence_No3d,Num_Sequence_No4d,Num_Sequence_No5d,Num_Sequence_No6d] = ...
	        ndgrid(NumRange_Gpoints_1d,NumRange_Gpoints_1d,NumRange_Gpoints_1d,...
			       NumRange_Gpoints_1d,NumRange_Gpoints_1d,NumRange_Gpoints_1d);% Note: Manual manipulation is related to Ndim
	Num_Gpoints_Nd = [ reshape(Num_Sequence_No1d,[length(NumRange_Gpoints_1d)^Ndim,1]) ...       % N dimensional Gauss Integration Points amount
					   reshape(Num_Sequence_No2d,[length(NumRange_Gpoints_1d)^Ndim,1]) ...       % Note: Manual manipulation is related to Ndim
			           reshape(Num_Sequence_No3d,[length(NumRange_Gpoints_1d)^Ndim,1]) ...
			           reshape(Num_Sequence_No4d,[length(NumRange_Gpoints_1d)^Ndim,1]) ...
			           reshape(Num_Sequence_No5d,[length(NumRange_Gpoints_1d)^Ndim,1]) ...
			           reshape(Num_Sequence_No6d,[length(NumRange_Gpoints_1d)^Ndim,1])];
	Num_Gpoints_Nd(find(sum(Num_Gpoints_Nd,2) < N_kLevel+1) ,:) = [];	        % N_kLevel+1 <= |Index_Set| <= N_kLevel+Ndim	   
	Num_Gpoints_Nd(find(sum(Num_Gpoints_Nd,2) > N_kLevel+Ndim) ,:) = [];
	Gpoints_Nd_SetHat = zeros(sum(prod(Num_Gpoints_Nd,2)),Ndim);                % the set of N dimensional Gauss integration points
% 	Gpoints_Nd_Set = zeros(sum(prod(Num_Gpoints_Nd,2)),Ndim);	                % the set of N dimensional Gauss integration points
	Weight_Nd_Set = zeros(sum(prod(Num_Gpoints_Nd,2)),1);                       % the set of N dimensional weight
	func_Value_Array = zeros(sum(prod(Num_Gpoints_Nd,2)),1);                    % Function value array corresponding to Gauss points
    [Deg_Sequence_No1d,Deg_Sequence_No2d,Deg_Sequence_No3d,Deg_Sequence_No4d,Deg_Sequence_No5d,Deg_Sequence_No6d] = ...
	   ndgrid(DegRange_1d,DegRange_1d,DegRange_1d,DegRange_1d,DegRange_1d,DegRange_1d);
	Degree_Nd_Set = [ reshape(Deg_Sequence_No1d,[(Polynomial_Order + 1)^Ndim,1]) ...
		           reshape(Deg_Sequence_No2d,[(Polynomial_Order + 1)^Ndim,1]) ...
				   reshape(Deg_Sequence_No3d,[(Polynomial_Order + 1)^Ndim,1]) ...
				   reshape(Deg_Sequence_No4d,[(Polynomial_Order + 1)^Ndim,1]) ...
				   reshape(Deg_Sequence_No5d,[(Polynomial_Order + 1)^Ndim,1]) ...
				   reshape(Deg_Sequence_No6d,[(Polynomial_Order + 1)^Ndim,1]) ];
    Degree_Nd_Set(find(sum(Degree_Nd_Set,2) > Polynomial_Order + Ndim),:) = [];
	Coefficient_Set = zeros(Num_monomials,1);
	counter = 0;
% The obtainment at N dimensional Gauss integration points
	for i_Row = 1 : size(Num_Gpoints_Nd,1)
		[Gpoints_Sequence_No1d,Gpoints_Sequence_No2d,Gpoints_Sequence_No3d,Gpoints_Sequence_No4d,Gpoints_Sequence_No5d,Gpoints_Sequence_No6d ] =...
			      ndgrid(1:Num_Gpoints_Nd(i_Row,1),1:Num_Gpoints_Nd(i_Row,2),...
			             1:Num_Gpoints_Nd(i_Row,3),1:Num_Gpoints_Nd(i_Row,4),...
						 1:Num_Gpoints_Nd(i_Row,5),1:Num_Gpoints_Nd(i_Row,6));            % Note: Manual manipulation is related to Ndim 
		Gpoints_Sequence_Nd = [reshape(Gpoints_Sequence_No1d,prod(Num_Gpoints_Nd(i_Row,:)),1) ...
			                   reshape(Gpoints_Sequence_No2d,prod(Num_Gpoints_Nd(i_Row,:)),1) ...
							   reshape(Gpoints_Sequence_No3d,prod(Num_Gpoints_Nd(i_Row,:)),1) ...
							   reshape(Gpoints_Sequence_No4d,prod(Num_Gpoints_Nd(i_Row,:)),1) ...
							   reshape(Gpoints_Sequence_No5d,prod(Num_Gpoints_Nd(i_Row,:)),1) ...
							   reshape(Gpoints_Sequence_No6d,prod(Num_Gpoints_Nd(i_Row,:)),1)];		
		Temp_Gpoints_Set = zeros(prod(Num_Gpoints_Nd(i_Row,:)),Ndim);
		for jndex = 1 : Ndim
			Gpoints_Noid = Chebyshev_Polynomial_Root(Num_Gpoints_Nd(i_Row,jndex));
			Temp_Gpoints_Set(:,jndex) = Gpoints_Noid(Gpoints_Sequence_Nd(:,jndex));
		end
		Gpoints_Nd_SetHat(counter+1:counter+size(Temp_Gpoints_Set,1),:) = Temp_Gpoints_Set;
		Temp_Sum_Degrees = sum(Num_Gpoints_Nd(i_Row,:));
		Weight_Nd_Set(counter+1:counter+size(Temp_Gpoints_Set,1),1) = (-1)^(N_kLevel+Ndim-Temp_Sum_Degrees)*...
		                                                            nchoosek(Ndim-1,N_kLevel+Ndim-Temp_Sum_Degrees)*...
																	1/size(Temp_Gpoints_Set,1);
		counter = counter + size(Temp_Gpoints_Set,1);
	end
% The Range [-1, 1] -> [a, b]: x = (a+b)/2+(b-a)/2*t
	Gpoints_Nd_Set = (a+b)/2 + (b-a)/2*Gpoints_Nd_SetHat;
% Function value obtainment at N dimensional Gauss integration points
	for i_Row = 1 : size(Gpoints_Nd_Set,1)
		func_Value_Array(i_Row,1) = ...
				 func(Gpoints_Nd_Set(i_Row,1),Gpoints_Nd_Set(i_Row,2),Gpoints_Nd_Set(i_Row,3),...
					  Gpoints_Nd_Set(i_Row,4),Gpoints_Nd_Set(i_Row,5),Gpoints_Nd_Set(i_Row,6));	
	end
% The coefficient obtainment
    % Root_Number = size(Gpoints_Nd_Set,1);
	for i_Order = 1 : length(Degree_Nd_Set)

			Temp_Vector = 2^(sign(Degree_Nd_Set(i_Order,1)-1)+sign(Degree_Nd_Set(i_Order,2)-1)+sign(Degree_Nd_Set(i_Order,3)-1)+...
			             sign(Degree_Nd_Set(i_Order,4)-1)+sign(Degree_Nd_Set(i_Order,5)-1)+sign(Degree_Nd_Set(i_Order,6)-1)).*Weight_Nd_Set(:,1).*func_Value_Array(:,1).*...
			     Numerical_Chebyshev_Polynomial(Degree_Nd_Set(i_Order,1)-1,Gpoints_Nd_SetHat(:,1)).*...
				 Numerical_Chebyshev_Polynomial(Degree_Nd_Set(i_Order,2)-1,Gpoints_Nd_SetHat(:,2)).*...
				 Numerical_Chebyshev_Polynomial(Degree_Nd_Set(i_Order,3)-1,Gpoints_Nd_SetHat(:,3)).*...
				 Numerical_Chebyshev_Polynomial(Degree_Nd_Set(i_Order,4)-1,Gpoints_Nd_SetHat(:,4)).*...
				 Numerical_Chebyshev_Polynomial(Degree_Nd_Set(i_Order,5)-1,Gpoints_Nd_SetHat(:,5)).*...
				 Numerical_Chebyshev_Polynomial(Degree_Nd_Set(i_Order,6)-1,Gpoints_Nd_SetHat(:,6));
			Coefficient_Set(i_Order,1) = sum(Temp_Vector);

	end

end